#include <buola/mat.h>
#include <buola/mat/prettyprint.h>
#include <buola/mat/markov.h>
#include <buola/mat/matrixops.h>
#include <buola/mat/statistics.h>
#include <buola/mat/random.h>
#include <buola/chart/cxychart.h>
#include <buola/chart/plots/clineplot.h>
#include <buola/chart/plots/cbarplot.h>

using namespace buola;

int main(int pNArg,char** pArgs)
{
    buola_init(pNArg,pArgs);
    
    random::engine().seed(0);
    
    const size_t lLen=1000000;
    const size_t lN=1000;

    //generate a random transition matrix
    mat::CMat_d lT=mat::random(lN,lN);
    msg_info() << sum(lT) << "\n";
    lT/=extend(sum(rows(lT)));
    msg_info() << sum(lT) << "\n";
    

    start_timer();
    mat::CVec_i lX=mat::generate_mm(lLen,lT);
    end_timer();
    
    msg_info() << mean(lX) << "\n";
   
    start_timer();
    mat::CMat_d lF=mat::estimate_mm(lX,lN);
    end_timer();
    
    //generate random transition matrix
    lT=mat::random(5,5);
    lT/=extend(sum(rows(lT)));
    mat::CRow_d lMu{0.3,0.2,0.1,0.22,0.18};

    //exercise 5
    {
        start_timer();
        while(true)
        {
            mat::CRow_d lSD=mpow(lT,64)(0,nAll);
            if(max(abs(lSD-lMu))<1e-8)
                break;
            
            for(int j=0;j<5;j++)
                lT(nAll,j)*=lMu(j)/lSD(j);
            lT/=extend(sum(rows(lT)));
        }
        end_timer();
        
        msg_info() << lT << "\n";
        msg_info() << mpow(lT,50) << "\n";
    }
    
    //exercise 8
    if(0)
    {
        const size_t lLen=10000000;
        mat::CVec_d lX=mat::zeros(lLen);
        
        for(int i=0;i<lLen;i++)
        {
            int lU=random::get<int>(0,2);
            
            if(lU==0)
                lX(i)=random::get<double>();
            else if(lU==1)
                lX(i)=pow(random::get<double>(),1.0/3.0);
            else
                lX(i)=pow(random::get<double>(),1.0/5.0);
        }
        
        mat::CVec_d lHCenters;
        mat::CVec_d lH=histogram(lX,20,lHCenters)*20.0/lLen;
        
        chart::PXYChart lChart=new chart::CXYChart;
        chart::PBarPlot lBar=new chart::CBarPlot;
        
        lBar->AddSeries(new chart::CSeries(lHCenters,lH));
        
        chart::PLinePlot lPlot2=new chart::CLinePlot;
        
        mat::CVec_d lPDFX=mat::linspace(50,0.0,1.0);
        mat::CVec_d lPDF=(1+3*pow(lPDFX,2)+5*pow(lPDFX,4))/3;
        lPlot2->SetSeries(new chart::CSeries(lPDFX,lPDF));
        
        lChart->AddPlot(lBar);
        lChart->AddPlot(lPlot2);
        lChart->Map(CSize_d(800,600));
        
        buola_main_loop();
    }
    
    if(0)
    {
        const int lN=10;
        mat::CVec_d lAlpha=mat::random(lN);
        lAlpha/=sum(lAlpha);
        const size_t lLen=1000000;
        mat::CVec_d lX=mat::zeros(lLen);

        for(int i=0;i<lLen;i++)
        {
            int lD=mat::sample(lAlpha)+1;
            lX(i)=pow(random::get<double>(),1.0/lD);
        }

        mat::CVec_d lHCenters;
        mat::CVec_d lH=histogram(lX,20,lHCenters)*20.0/lLen;
        
        chart::PXYChart lChart=new chart::CXYChart;
        chart::PBarPlot lBar=new chart::CBarPlot;
        
        lBar->AddSeries(new chart::CSeries(lHCenters,lH));
        
        chart::PLinePlot lPlot2=new chart::CLinePlot;
        
        mat::CVec_d lPDFX=mat::linspace(50,0.0,1.0);
        mat::CVec_d lPDF=0*lPDFX;
        for(int i=0;i<lN;i++)
            lPDF+=lAlpha(i)*(i+1)*pow(lPDFX,i);
        lPlot2->SetSeries(new chart::CSeries(lPDFX,lPDF));
        
        lChart->AddPlot(lBar);
        lChart->AddPlot(lPlot2);
        lChart->Map(CSize_d(800,600));
        
        buola_main_loop();
    }

    {
        const size_t lLen=1000000;
        mat::CVec_d lX=mat::zeros(lLen);
        
        for(int i=0;i<lLen;i++)
        {
            double lY=-log(1.0-random::get<double>());
            lX(i)=pow(random::get<double>(),1.0/lY);
        }

        mat::CVec_d lHCenters;
        mat::CVec_d lH=histogram(lX,20,lHCenters)*20.0/lLen;
        
        chart::PXYChart lChart=new chart::CXYChart;
        chart::PBarPlot lBar=new chart::CBarPlot;
        
        lBar->AddSeries(new chart::CSeries(lHCenters,lH));
        
        lChart->AddPlot(lBar);
        lChart->Map(CSize_d(800,600));
        
        buola_main_loop();
    }
    
    return buola_finish();
}
